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SUMMARY 


A hybridized method for computing the flow about blunt bodies is presented. 
In this method, the flow field is split into its viscid and inviscid parts. 

The forebody flow field about a parabolic body is computed. For the vis- 
cous solution, the Navier-Stokes equations are solved on orthogonal parabolic 
coordinates using explicit finite differencing. The inviscid flow is determined 
by using a "Moretti type" scheme in which the Euler equations are solved, using 
explicit finite differences, on a nonorthogonal coordinate system which uses the 
bow shock as an outer boundary. The two solutions are coupled along a common 
data line and are marched together in time until a converged solution is 
obtained . 

Computed results, when compared with experimental and analytical results, 
indicate the method works well over a wide range of Reynolds numbers and Mach 
numbers . 


INTRODUCTION 

Over a number of years, highly refined techniques have been developed for 
obtaining flow- field solutions about blunt entry bodies at supersonic and hyper- 
sonic speeds. Early techniques included the inviscid shock-layer method of 
Moretti and Salas (ref. 1) and the viscous shock-layer method of Davis (ref. 2) 
and others. The more recent work of Kumar (ref. 3) has extended viscous blunt- 
body calculations to include solutions for small angles of attack and for the 
plane of symmetry. In one way or another, all of these methods use a coordinate 
system tied to the body shape along with a transformation of the physical space 
into a uniform rectangular mesh in the computational plane. The dependence of 
the transformations on body shape necessarily limits these techniques to flow 
over smooth shapes, since any attempt to carry the solutions past a sharp 
forebody/afterbody juncture will result in a singularity in the coordinate 
transformation . 

The ability to compute the complete flow field about a finite entry body 
has always been a desirable but, until recently, unrealized goal due largely to 
inadequate computer resources. Recently, Peyret and Viviand (ref. 4) and Gnoffo 
(ref. 5) have solved the complete flow about finite analytical bodies at an 
angle of attack of zero degrees. In both papers, Navier-Stokes solutions were 
obtained on a physical mesh which is tied to the body shape. The solution of 
Peyret and Viviand is for a body having a sharp forebody/afterbody juncture, 
while Gnoffo uses a smooth, analytical body shape. Both methods also use shock 
capturing to handle the bow shock. The shock capturing technique works well for 
flows at low Reynolds numbers. However, at high Reynolds numbers, the number of 
mesh points required to adequately resolve the shock are prohibitive. Gnoffo, 
in unpublished work at NASA Langley Research Center (LaRC) , has floated a 



discrete shock in his mesh to handle these conditions at high Reynolds numbers 
while solving the full Navier-Stokes equations everywhere in the flow field. 

In this paper, we will demonstrate a technique which simultaneously solves 
for the flow in the inner, viscid region and the outer, inviscid region about 
the forebody of an axisymmetric parabola at 0° angle of attack. This tech- 
nique, which is directly applicable to the finite body problem, will minimize 
the number of grid points required in the computational domain. Solutions for 
the outer, inviscid region are found using a "Moretti type" method which is 
tied to an arbitrary line through the flow field rather than to the body sur- 
face. The inner, viscous solution is found by solving the full Navier-Stokes 
equations over a region bounded by the body surface and the lower bound of the 
inviscid solution. The solution technique has been structured for the Control 
Data STAR-100 computer at Langley Research Center (ref. 6) , thus, the method 
will be explicit in nature. 


SYMBOLS 


v 


/\ j/\ 

local speed of sound, a /U 
constant in equation (2) 
specific heat at constant pressure 
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(eq. (A3)) 


nose radius 
body surface length 

AA /A 

time, tU /L 

’ COf 

/\ J/\ 

temperature, T/T 
velocity, 
x-component of velocity 
p-component of velocity 
velocity, v/ (eq. (A3)) 

y-component of velocity 
^-component of velocity 
Cartesian coordinate direction 
Cartesian coordinate direction 
ratio of specific heats 
increment of length 


^max % 


parabolic coordinate 
computational step size, At/A 
viscosity, 

parabolic coordinate 
density, p/p m 
(J) transformed coordinate 

Subscripts : 


common 


normal to shock 


I 


o 


reference to body surface 


s reference to post- shock conditions 


t 

= 3( 

)/3t 

w 

wall 


n 

= 3( 

)/3n 

S 

= 3( 

)/3£ 


00 reference to free-stream conditions 

Other notations : 

Arrow over symbol (e.g., F) denotes vector. 

Bar over symbol (e.g., t) denotes coordinate in inviscid flow field. 
Circumflex over symbol (e.g., u) denotes dimensional variable. 


ANALYSIS 


Several factors were taken into account while determining the appropriate 
coordinate system to be used with the computational technique described in this 
paper. The shock capturing technique was not considered because of the large 
number of grid points required to adequately resolve the bow shock at the high 
Reynolds numbers of interest. Also, the shock floating technique was not con- 
sidered because it requires a surplus of mesh points in the computational field 
as illustrated in the following sketch and because the logic required to float 
the bow shock in the grid would nullify the computational advantages of vector 
processing. Thus, it was determined that the physical plane to be covered 
should be an envelope bounded by the bow shock, the body, the outflow boundary, 
and the stagnation line as shown: 



OUTFLOW 

BOUNDARY 


Rectangular Computational Plane 
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The global coordinate system consists of two distinct but interdependent 
local coordinate systems. Both systems are constructed so that, through a 
coordinate transformation, the numerical computations are carried out on a uni- 
form rectangular mesh. 


Local Coordinate System for the Viscous Region 

The physical coordinate system used in the viscous computations is the 
orthogonal parabolic system used by Peyret and Viviand (ref. 4) and is formed by 
the intersection of curves belonging to two cofocal parabolas as shown in fig- 
ure 1(a). The parametric equations for the Cartesian coordinates in this system 
are 



As seen in figure 1 (b) , the axisymmetric body surface is formed by rotating one 
of the coordinate lines, specified by a constant value of p = p , about the 
x-axis. 

The computational technique, to be described in the section "Viscous 
Solution," requires a mesh ratio of one in the computational plane. Shown in 
figure 1(c) is the physical plane transformed to the uniform computational plane 
which is divided into equal increments of E and p (i.e., A£ = Ap) . The 
mesh in the physical plane expands naturally in both the E,~ and p-directions 
and provides the greatest spatial resolution in the region of the nose in both 
the E~ and p-directions and along the wall in the p-direction. Since the size 
of the viscous region for any given conditions is fixed, the spatial resolution 
can be controlled by either increasing or decreasing the incremental size of E, 
and p which either decreases or increases the number of grid points. Con- 
trolling the resolution in this manner is sufficient when working at low 
Reynolds numbers, but it is inadequate at high Reynolds numbers because of the 
large number of grid points needed when very small values of Ap are needed 
near the wall. The problem was overcome by applying a Blottner (ref. 7) type 
transformation to the p-coordinate in the form 


I a X/AX _ x \ 

L^x . J 


+ n 


o 


( 2 ) 


Now, computations are performed in the XrE grid where AX = A^ . This allows 
reasonable values of A E, to be chosen while the proper distribution of points 
in the p-direction is controlled by X max r Ap, and the choice of the con- 
stant A. The stretched £,p grid is shown in figure 1(d). 
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Local Coordinate System for Inviscid Region 


The coordinate system used for the inviscid computation is similar to that 
used by Moretti and Salas (ref. 1) except that here the lower boundary of the 
grid system is tied to a line of constant r] (the common line p c in 
fig. 1(b)) rather than to the body surface. A new set of transformed coordi- 
nates (5, T"Lt) is defined by the transformation equations 


5 = 5 


n - n r 


n«.(5,t) - n. 


t = t 


(3) 


which transform the region between the common line r) c and the shock wave, 
n s = r) s (£,t), into a rectangular region with the shock forming the outer bound- 
ary. (See fig. 2.) The grid spacing in the ^-direction and the time t are the 
same for the inviscid and viscous regions. However, the coordinate system in 
the inviscid region is nonorthogonal because of the transformation used for r) , 
and since the shock location varies with time (until steady state is reached) , 
the coordinate system also varies with time. 


The Global Coordinate System 

As illustrated in figure 3, the two meshes in the physical plane are joined 
along a constant grid line, r) = f| . Information along this line is computed by 
the viscous portion of the code and serves as a boundary condition for the 
inviscid portion of the code. Data on the outer boundary of the viscous solu- 
tion, g = ^interpolated' are ^ oun< ^ ky interpolating along lines of constant E, 
in the inviscid solution. The way in which the flow field is divided into an 
inner and outer solution will be discussed in a later section. 


Computations are carried out using velocities in the parabolic rather 
than Cartesian system as shown in figure 4. The physical velocities along 
E, - constant and r| = constant are as follows: 


U 


-nu + £v 

h 


(4) 


{/ 


E, u + gv 

h 


(5) 


Now, set h = h^ = h 2 


(e 2 4. n 2 ) 1/2 


and define a set of coinvariant velocities 


u = hU and v = hV 


(6) 
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The use of these coinvariant velocities in the governing equations eliminates 
errors caused by the differencing of nonlinear metric coefficients. 


Viscous Solution 

The solution in the viscous portion of the flow is obtained by integrating 
the full, compressible, time-dependent Navier-Stokes equations on the orthogonal 
parabolic grid described in the previous section. The integration is carried 
out by applying a modified version of the two-step Richtmeyer differencing 
algorithm described by Thommen in reference 8 to the appropriate system of 
equations . 

The governing equations are written in the conservative form 


-> 



9f (W) 

3 ? 


3g(w) 

9n 


+ H(W) 


+ S(W) 


(7) 


where the terms in equation (7) are the conserved quantities W, the convective 
terms F and G, the geometry terms H arising from a non-Cartesian coordinate 
system, and the viscous dissipation and heat conduction terms S. The specific 
form of the Navier-Stokes equations used for this two-dimensional, axisymmetric 
flow in an orthogonal parabolic coordinate system is developed in appendix A. 

The difference algorithm in equations (8) and (9) is used to integrate 
equation (7) on a nine point computational module for any W. The variable W 
is defined in the time plane t at the grid points represented by circles and 
in the plane t + At/2 at the points marked by crosses in this sketch: 
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The two-step differencing scheme is written as follows: 


->t+At/2 
W ±l/2 , 0 


1 Aft -ft \ At f/f- \t /•% 

2 \ ± 1 , 0 Vo) 2 |{ F ?:)±l/ 2 / 0 + ( G l 

1 Ait :+t \ 1 //t -+t \ 

2V ±1,0 0,o) 2 \ ±1 , 0 0,o) 


\ fc 

n/±i/2,o 


(8a) 


w 


:t+At/2 
0 , ± 1/2 


2 \ 0 , 


->-t 

, + w 

±1 0,0 


At 

2 


WU * K 


<Vo,±l/2 + <S)o.±l/2 

.0)- V: 


.±1 + *o,o 


(8b) 


-*t+At -*t . 

w „ „ = w„ „ - At 

0,0 0,0 


«* C >oIf 72 * <%)olo t/2 - * 


H o,o s o 


, 0 ) 


(9) 


The first-order functions in equations (8) and (9) can be written for 
either F or G as 



t 

±1/2,0 







<0: 


& 


n/±i/ 2 ,o 4AnV o,+i 


+ *±i,+i 


- ? 
Or 


- F fc 

± 1,-1 


( 10 ) 


/t, yt+At/2 = 1 A/t+At/2 :tt+At/2\ 

\ C/0,0 AC \ +1/2,0 “ -1/2, 0/ J 

Their counterparts (g^^, <W ±1/2 , and are found by 

permuting indices in equations (10) . The second derivatives used in the viscous 
term S are defined below for the general variable ip: 
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<K>5>o,o 



b 


t 

+ 1 / 2,0 




+ 1,0 



b- l/2 , 0 




( (b ^5 ) n)o,0 4A?Ari b o, +1/2 (^+1,0 + t-1,+1 ^-1,0 ^-l,+l) 


- b 


t 

0 , - 1/2 




+ 1,0 







( 11 ) 


b 


t 

± 1 / 2,0 



+ b 


± 1,0 




Thommen first recommended that the H and S vectors be reevaluated at 
t + At/2. He later found this was unnecessary for Because of the large 

variation of the metric coefficients near the stagnation line for this particu- 
lar coordinate system, we have obtained better results in the stagnation region 
by not reevaluating H. 


Inviscid Solution 

The solution in the inviscid portion of the flow is obtained by integrating 
the compressible, time-dependent Euler equations in the transformed coordinate 
system described in a preceding section. The integration is carried out using 
the two-step MacCormack finite-difference algorithm (ref. 9) which is modified 
by using one-sided differences at the shock boundary. 

The conservative, vector form of the Euler equations written in orthogonal 
parabolic coordinates (5,n»t) is 


-> 



= -Fr (W) - G n (W) + H(W) 


( 12 ) 


where W, F, G, and H are the same vector quantities defined in appendix A. 
By using the transformation equations (eqs. (3)), the system of equations 
(eq. (12)) can be transformed from parabolic coordinates (5,1"), t) to a set of 
coordinates in the computational space (5,r|,t). Thus, equation (12) becomes 


W- 

t 


9f) ■+ 
— 1 w- 
9t n 


-* 



In 

95 n 


3n r _ 
3n n 


-> 
+ H 


(13) 


Since only the steady-state solution is of interest, the term ttt W- is 

d t ri 

neglected and equation (13) becomes 
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(14) 


-> 




3£ n 


3 q r._ 
3n n 


+ H 


The details of this transformation are given in appendix B. 

The difference algorithm in equations (15) and (16) is used to integrate 
equation (14) on the computational module illustrated below: 



The variable W is defined in the plane t = 0 at the points represented by 
the circles in the sketch above. The two-step MacCormack difference scheme is 
then written as follows for the solution of w at the time plane t + At: 

Predictor step - 


->t+At -n. 

W = - At 

0,0 0,0 


3F\ fc + /3n\ fc / 3 f f 

3 ?/o,o v?/o,o\3n/ f 


0,0 


rtf /9G Nt 

9n ^o,o\3n 


0,0 


-y t 

- (H) o,o 


(15) 


Corrector step - 


->t+At 1 7-H; ->t+At 

W = — < W + w -At 

0,0 2 \ 0,0 0,0 


3f 

35 


t+At 


rt) t+At hl) t+lst 

35/0,0 \3o/o,o 


_ (S) 5Stl 

\3n/o,o \3fi/o,o °'° 


(16) 
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At the interior points on the computational grid, two-point backward differ- 
ences are used to replace the derivatives of F and G m the predictor step 
and two-point forward differences are used in the corrector step. For the grid 
points along the bow. shock, three-point backward differences are used in both 
the predictor and corrector steps. The procedure used to track the movement of 
the bow shock wave will be described in the next section. 


Solution for Shock Wave 

A method similar to that used by Tannehill, Holst, and Rakich (ref. 10) has 
been used to track the movement of the shock wave during the transient portion 
of the solution. The method can be summarized as follows (details are given in 
appendix C) . First, the pressure is computed on the downstream side of the 
shock wave using the two-step difference scheme described in the preceding sec- 
tion with three-point backward differences used to compute 3 f/3r and dG/dri 
at the shock wave. With the pressure at the shock wave P s known, the density 
is computed using the equation 


(17) 




1 

Ps 


+ f (P s 




'2 \ _ l (p _ 

2 s 


Poo) 


Then, the shock velocity is computed using the equation 


V 


s 


(Voo 



(18) 


where V~ and (V„) are defined in sketches (c) and (b) in appendix C. 

° 11 OO 

Next, the movement of the shock wave (see fig. 5) is computed using the follow- 
ing two-step difference method: 


Predictor step - 


, . t+At , , t . 

(n s ) = (n s ) - At 


V, 


h cos t, c 


Corrector step 


(n„) 


t+At 


, , t At 

= (n«) - 


' V s 
ih cos 


, t+At 


h cos 


(19) 


( 20 ) 
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where 


= 


tan 



(21) 


The derivative 3r) s /3£ in equation (21) is approximated by central differences 
in both the predictor and corrector steps. After the new location of the shock 
wave has been determined, the properties (except pressure) on the downstream 
side of the shock are computed using a transient shock-wave analysis. 

As the solution approaches convergence, the shock velocity approaches zero 
and the shock wave remains essentially fixed in space. This usually occurs 
before the remaining properties in the flow field converge. The shock veloci- 
ties at convergence are generally less than 10“^ in magnitude. 


BOUNDARY AND INITIAL CONDITIONS 


Boundary Conditions 


Viscous regio n.- Along the body surface, the usual no-slip boundary condi- 
tions, u = 0 and v = 0, are imposed. Wall pressure is found by setting the 
normal pressure gradient 3P/3r| equal to zero while the wall temperature is 
held constant at some predetermined value. Thus, the wall density is found 
from the equation of state. 

Symmetry conditions 


p(n,-A£) = 
u(u,-A£) = 
v(r|,-A£) = 
e(u,-AC) = 
P(n,-A£) = 


P(U,A£) 
u(n,A£) 
-v(ri,A£) \ 

e (n,A£) 
P(n,A£) 


( 22 ) 


are imposed along the stagnation line E, = 0. At the outflow boundary, a 
second-order extrapolation of the conserved variables 


W(£ 


'max ,r ^ 3W ^max A£,ri) 3w ^ max 


2A£,n> + w (5 max " 3A£,U) 


(23) 


is used. 


12 



As shown in figure 3, the inflow boundary for the viscous solution is 
determined by linearly interpolating along lines of constant E, on information 

in the inviscid solution to get values along n = ri . . . , _ . 

3 'interpolated 

In viscid solution .- As shown in figure 3, the lower boundary of the invis- 
cid solution is fixed by making it coincident with the line p = n c , the last 
line of computed information in the viscous solution. Along the stagnation 
line, the symmetry conditions are the same as those used in the viscous 
solution . 

Values along the outflow boundary are determined by a second-order extrapo- 
lation of values at the interior points, while the values on the inflow boundary 
are determined by the unsteady shock relations discussed in the inviscid flow 
section . 

Initial conditions .- Initial conditions are imposed in the following 
manner. The wall surface temperature is fixed and velocities at the wall are 
set equal to zero while a Newtonian pressure distribution is determined for the 
given body shape. An initial shock shape, which is usually a line of con- 
stant ri , is picked. Steady shock relations are then used to determine flow 
conditions behind the shock. The pressure and velocities are then linearly 
interpolated between the shock and body along lines of constant E,. Finally, 
the temperature along each line of constant E, is set equal to the shock tem- 
perature up to, but not including, the wall. The density distribution is then 
found using the equation of state. 


Numerical Stability 

The inviscid numerical stability limit for the two methods being used is 
the well-known CFL condition 


At 


c o/c 

M. , M . 1 1 j_ 

Ax Ay a VAx Ay 


(24) 


where C' =1 for the MacCormack scheme and C' =2 for the Richtmeyer scheme. 
For both methods, the Courant number Cq is taken to be 0.95. When equa- 
tion (24) is applied to the present coordinate system and u and v are used 
as defined in appendix A, the equation becomes 


At 

A 


(C 0 /C')h 2 
| u [ + | v + \f 2 ha 


(25) 
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Global Operational Characteristics 

The preceding procedure has been programed and run on the Control Data 
STAR-100 computer at Langley Research Center. The computational mesh size used 
depends on the flow Reynolds number and on the need to adjust the relative posi- 
tion of points in the two meshes. Inviscid mesh size varied from 6 x 21 to 
12 x 21, and the viscous mesh varied from 14 x 21 to 25 x 21. Normal mesh 
spacing in the viscous region was determined by requiring that the stagnation- 
point cell Reynolds number equal two. 

Initially, solutions were obtained using a single value of At which was 
determined by finding the minimum of all local values of At. However, the 
combination of a high Reynolds number with the restriction on cell Reynolds 
number led to extremely small values of At near the wall. This led us to try 
using the local minimum At at each point, as reported by Kumar (ref. 3). This 
procedure yielded good results and gave us a reduction of approximately 50 per- 
cent in the number of iterations required to reach convergence for the highest 
Reynolds number runs. 

The actual running of a computer code as described always raises a number 
of operational problems which cannot be quantified. In particular, a number of 
problems have appeared at the juncture of the inviscid and viscid flow meshes. 
Each problem has been resolved, but for most, there is a lack of a definitive 
cause-and-ef f ect relationship. 

The first and most obvious problem is determining whether the flow field is 
properly split into its viscous and inviscid components. We quickly found that 
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the code simply would not run if the inviscid region extended into the viscous 
region. We surmounted this problem by being very conservative in our estimate 
of the extent of the viscous region. 

We also found that, in order to maintain a stable computation, the inter- 
polated row of data in the viscous solution must be coincident with, or have a 
value of r| greater than, that of the second row of data in the inviscid 
solution . 

The flow-field solution was considered to be converged when the relative 
convergence criterion 


t+At t 
P ~ P 


^ 10 -5 


was satisfied at every point in the flow field. 


RESULTS AND DISCUSSION 

The preferred method of validating numerical techniques is comparison with 
experimental data. Experimental data for flow over parabolic bodies are scarce 
and basically limited to drag coefficients. We have used the data of Little 
(ref. 11) and information generated in the Langley continuous-flow hypersonic 
tunnel (CFHT) for drag data. Shock shapes were obtained from tests run in the 
Mach 6 high Reynolds number tunnel at the Langley Research Center. Also, com- 
parisons are made with established numerical methods to determine the accuracy 
of pressure and heat-transfer distributions. All computations have been made 
using the Sutherland viscosity formula and a perfect gas with a ratio of 
specific heats of 1.4. 

Our first concern was determining the effect of splitting the flow field on 
the overall accuracy of the solution (specifically, determine whether the physi- 
cal joining of two different computational meshes would distort the solution in 
the vicinity of the juncture) . To check this, we compared the complete inviscid 
solution with the coupled solution for two different sets of free-stream flow 
conditions. For the first condition, shown in figure 6, N R _ = 2 x 10® and 

■^oo 

= 5.87. Here, temperature profiles are plotted for the stagnation line 
(£ = 0) and a point approximately halfway around the body. Both curves indicate 
a smooth transition from one mesh to the other at the point marked "mesh 
juncture" which is the last computed point in the viscous solution. Also, the 
two solutions agree very well down to the point where boundary-layer effects 
dominate the flow. However, there is, as expected, some difference between the 
two curves due to the effect of boundary-layer growth on effective body shape. 
The plot presented in figure 7 is similar except that the N R „ has been 

reduced to 3 x 10® and = 9.92. Again, good results were obtained, but the 

difference between the two curves at the downstream location is more pronounced 
due to the more rapid boundary-layer growth. 
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Next, a comparison is made with experimental drag coefficients. In fig- 
ure 8, numerical results are compared with the data of Little (ref. 11). The 
present method can be used only for body configurations which have supersonic 
outflow boundary conditions outside of the boundary layer. There is, however, 
good agreement with experimental data over the range of conditions which can be 
duplicated by the present program. It is also interesting to note the viscous 
effect on the drag coefficient, even at these relatively large Reynolds numbers, 
as indicated by the divergence of the viscid and inviscid curves as body length 
increases . 

The second comparison is made with experimental drag coefficients based on 
data taken in the Langley CFHT and is shown in figure 9. Here, the comparison 
of variations in the drag coefficient is made for a range of Reynolds numbers 
on a 30° paraboloid. 

Our last comparison is made with experimental data for the shock shape 
shown in figure 10. The experimental shock shape was obtained from schlieren 
photographs taken in the Mach 6 high Reynolds number tunnel at the Langley 
Research Center for the conditions M^, = 5.8725 and N Re = 2 x 10^. Again, 

there is excellent agreement between experimental and numerical results. 

There is no experimental heat- transfer data available for comparison with 
results generated by our present method. Therefore, we have chosen to make our 
comparisons with the results generated by several proven computer codes. For 
M„ = 10.34 and N R _ = 5 x 10^, results are shown in figure 11. The results 

from the present method are compared with computations using the Navier-Stokes 
solution of Graves (ref. 12), the viscous shock-layer solution of Kumar 
(ref. 3) , and the boundary- layer solution of Anderson (ref. 13) which uses the 
inviscid edge conditions from Marconi (ref. 14) . As can be seen in figure 11, 
there is extremely good agreement among all methods. In figure 12, similar 
results are shown for the free-stream conditions M^ = 5.8725 and 
N R = 2 x lo 6 . 

Kc co 


CONCLUDING REMARKS 

In the previous sections, we have investigated the viability of a hybrid 
viscous/inviscid solution technique for determining the flow about an axisym- 
metric parabola. In this investigation, we have found that it is possible for 
a flow-field computation to be marched in time to convergence, where the field 
is split into an inviscid region governed by the Euler equations and a viscous 
region governed by the Navier-Stokes equations. Furthermore, we have experi- 
enced no problems at the juncture of the two flow regimes when the flow field 
has been properly split into its viscous and inviscid parts. By comparison 
with both analytical and experimental results, we have demonstrated that the 
technique works well over a wide range of Reynolds numbers and Mach numbers. 
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As presented, this technique reduces the number of grid points required 
in the computational field, and by its explicit nature, the method is well 
suited to take advantage of the high computational speed of vector processing 
computers . 

Although presented for forebody flows only, there are no conceptual prob- 
lems to extending this technique to the flow about finite, sharp-cornered 
bodies. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
July 6, 1979 
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AXISYMMETRIC NAVIER-STOKES EQUATIONS IN 
ORTHOGONAL PARABOLIC COORDINATES 
General Equations 

The dimensional Navier-Stokes equations , ^ written in conservation form for 
a general, curvilinear coordinate system, may be stated in terms of the system 
metric coefficients, h^, h 2 , and h 3 , as follows: 

Continuity - 

3 8 8 3 

g^-(h 1 h 2 h 3 p) + 3^'(h2 h 3 pu l ) + g^' (h l h 3P u 2 ) + g^j(h 1 h 2 pu 3 ) = 0 (Ala) 


x-, -momentum - 


3t (h l h 2 h 3 pu l ) + 3 Xl 


h 2 h 3 (p + p U 2) + ^-(h^pu^) + g~-(h 1 h 2 pu 1 u 3 ) 


/ 2 \ 9h 2 3h^ 9h^ / 2 \ 9h 3 

- ^p + pu 2 Jh 3 + P u ! u 2 h 3 3 ^ + P u l u 3 h 2 9^ - [V + P u 3 ) h 2 9^ 


“ 9~ (h 2 h 3 T ll ) + 9x 2 (h l h 3 T 21 ) + 9x 3 (h l h 2 T 31 ) 

9h]_ 3hi 9h2 9h3 

+ x 21 h 3 9 ^ + T 31 h 2 9 ^J ' x 22 h 3 9 ^ " X 33 h 2 9 ^ 


(Alb) 


^For convenience, the superscript denoting dimensionality ~ has been 
omitted from these equations. 
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x 2 -momentum - 


9 9 9 ” / 2 

9^(h 1 h 2 h 3 Pu 2 ) + g^T'( h 2 h 3 pu l u 2 ) + 777 _ h l h 3 ( pu 2 + P )_ 


pu 2 + P )| + 3x 3 (h l h 2 pu 2 u 3 ) 


/ 2 \ 9 h 2 9 h 2 / 2 \ ^^1 

- (PU3 + P)h X g^~ + PU 2 U 3 h X g^~ + PU lU2 h 3 g^~ - (PU X + P)h 3 


9 9 9 

37^' (h 2 h 3 T 12 ) + 37^ (h 3 h l T 22 ) + '3^' (h l h 2 T 32 ) 


3 h 2 9 h 2 9 h 3 3 h-^ 

+ T 23 h l g7~ + x 12 h 3 9^ ‘ X 33 h l 97^ " X H h 3 9^ 


x 3 -momentum - 


|^(h 1 h 2 h 3 pu 3 ) + ~(h 2 h 3 p Ul u 3 ) + g|-(h 1 h 3 pu 2 u 3 ) + 4[ hl h 2 (pu^ + pj] 


12 \ ^hi 0113 0113 / 2 \ 

- (p U;L + pjh 2 g— + pu l u 3^2 377 + pu 2 u 3 h l 977 - \ Pu 2 + p j h l 


9h 


9h- 


9h- 

9x- 


9 9 9 

97^ (h 2 h 3 x 13 ) + 97^ (h 3 h l X 23 ) + 7x^ (h l h 2 T 33 3 


3h 3 3h3 3h-^ 3h 2 

+ X 31 h 2 97^ + T 23 h l 97^ “ X ll h 2 97^ ~ X 22 h l gTJ 


Energy equation 


9t (h l h 2 h 3 E) + 9 X;L [ h 2 h 3 u l (E + P) J + 3x 2 &l h 3 u 2 (E + P G + 9x 3 [ h l h 2 u 3 (E 


3 / h 2 h 3 3t \ 3 / h l h 3 3t \ 3 / h l h 2 „ 3t 


3x^1 h-^ K 3 x-ji ~ r 3 x 2 l h 2 ^ 3x 2 ) T 3x 3 ^ h 3 ^ 3x 3 


K 


K 


3x- 


-[ h 2 h 3 (u 


l x ll + U 2 X 12 + U 3 X 13 ) ] + 3x- 


-[hih 3 (u 


l x 12 + U 2 X 22 + u 3 


+ 37r[ h l h 2 (u l x 13 + U 2 X 23 + U 3 X 33 ) ] 


(Ale) 


(Aid) 


+ P)] 


X 23 ) ] 

(Ale) 
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where E is defined as p^e + 
defined as follows: 


u x + vi 2 + u 3 


The viscous stress terms are 


T 11 = 


(l 8u l U 3 8h l 

9x-^ + h^h 2 9 x 2 


where 


XV • V = 


2y 


T 22 2y lh 2 3 x 2 


3h l h 2 h 3 


1 ^ u 2 . u 3 3h 2 


t 33 - 2 V\u_ 


1 8u 3 


t 23 - V 


13 


I 3 -A 

h 2 3x 2 V h 3 i 


12 



h-^ 9x^ 

h 

^2 9 | 

/u 2 ’ 

h^ SXgl 

V h 2; 

h 3 9 | 

^ u 3 

h l 3x l 

\ h 3, 

^1 _9_| 

f u l 

h 2 3 x 2 ! 

l h l 


h 3 h l 

9x 3 

l) + 

3 

3x 2 l 

U 1 

9h x ' 

h l h 2 

9x 9 

u 2 

9h 3 ' 


+ XV • v 


+ XV • V 


+ XV • V 


t 21 “ t 12 


t 31 T 13 


t 32 “ t 23 


Equations in Parabolic Coordinates 


For the orthogonal parabolic coordinate system, the following substitutions 
are made in equations (Al) : 
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X 1 = 5 

x 2 

= n 

x 3 

= <f> 

uj = V 

u 2 

= u 

u 3 

= W 


hx = (? 2 + n 2 ) 1/2 h 2 = a 2 + T] 2 ) 1/2 h 3 = 


Furthermore, W and 9/9(j) are set equal to zero for axisymmetric flow. These 
substitutions and restrictions, when applied to the general equations, lead to 
the following set of equations : 

Continuity equation - 


9 t (hih 2 h 3 p) + gg.(h 2 h 3 pl/) + 3^ (h 1 h 3 p(i) 0 


^-momentum equation - 


(A2a) 


^(hxh 2 h 3 pl/) + |^[h 2 h 3 (pt/ 2 + P)] + |^(hxh 3 pUl/) 

9h_ 9h. 9h_, 

- (pU 2 + P) h 3 _ + p ui/h 3 - Ph 2 

3 3 

= 9^' (h 2 h 3 T ll ) + 97i' (h l h 3 T 21 ) 

9h^ 9h 2 9h 3 

+ X 21 h 3 9r| T 22 h 3 9| X 33 h 2 ' 

T|-momentum equation - 

|^(hxh 2 h 3 pU) + ~(h 2 h 3 pUt/) + ^[hxh 3 (pu2 + P)] 

3h^> 3hp 0 3hi 

Ph i aTj - + P ul/h 3 g| (p 17 + p)h 3 aTp 

3 3 

= 3|-(h 2 h 3 x 12 ) + 3^(h 3 h 1 T 22 ) 

9h 2 9h 3 9hx 

+ X 12 h 3 9| X 33 h l 9ri x ll h 3 9VP 


( A2b) 


(A2c) 
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Energy equation - 

fe(h 1 h 2 h 3 E) + ^[h 2 h 3 l/(E + P)] + |^[h 1 h 3 U(E + P)] 


3 ( h 2 h 3 

3t\ 

3 / h l h 3 

3t' 

1 


3q h-L 

h) 

+ 9 n^ h 2 

K 

3ri j 

1 


3 r 

+ 8^L h 2 h 3 

1 Vt 11 

+ Ut 12 )] 

3 1 

+ an 1 

h l h 3 ^ T 12 + ^ T 22^j 

(A2d) 


Stress terms - 


T 11 = 


and 


XV • V = 


31/ 

U 

3£ h 1 h 2 

2y 

'3 

3h l h 2 h 3 

L 9 ^ 




+ XV • V 


-(h 2 h 3 l/) + g^(h 1 h 3 (J) 




( 1 dU 1/ 3h x 

t 22 ~ li- " + 


t 33 - 2y 


h 2 3n h i h 2 3? 

1/ 3h 3 tj 3h 3 

h 3 h^ 3? h 2 h 3 9r l 


+ XV • V 


+ XV * V 


t 23 ~ 0 


t 13 - 0 


: 12 - ^ 


h2 9_/_U_\ ^ L-(JL 


L h l 9 ^\ h 2/ h 2 9r A h l/J 


J 


(A2e) 


When the coinvariant velocities 


u = hU 

and 

v = hV 


(A3) 


are introduced into the equations (A2) , the result is as follows: 
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Continuity equation - 


h lh2p) = ~ {pv) 5 " (pu) n ' p (^ + f) 


5-momentum equation - 


-j^-(h 2 pv) = - -^(pv 2 + h 2 P) - -^pj-(puv) + Jj>h 2 P + p (u 2 + v 2 )j 

/u v \ h[/3 3 \ 

' pv (n + c) + T u + ST T 2ij 

/ 5 2 + h 2 \ . (h 2 + 2r| 2 \_ 5 h 

1 Ch ) T 11 + I nh ] T 21 “ h T 22 - £ t 33 




p-momentum equation 


g^(h 2 pu) 


|^-(puv) - |^-(pu 2 + h 2 P) + j^2h 2 P + p (u 2 + v 2 )J-^- 


■Cn 


v \ h 

- pu( - + ^ + - 


3t 12 3 t 22' 
35 + 3n ) 


2l 2 + h* . 

+ 'Eh l T 12 


nh 


h n 

22 - n T 33 h T 11 


Energy equation - 


3t 


(h2E) = - fe[ v(E + p) ] - U u(E + p 3 - (e + p) (^ + f) 


+ U K ^ 

35 1 35 


\ 3 / 3t\ fl 3t l 3 t\ 1 , \ 

) + 3n\ 3n j + K (j 35 + n 3n j + 5 (vT n + ut i 2) 


+ n (vT 12 + uT 22 ) + 3^ (vT 11 + uT 12 ) + 3r| (vT 12 + uT 22 ) 


(A4a) 


(A4b) 


(A4c) 


(A4d) 
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Stress terms - 


and 


T n = 7l( h2 || - vC + unj + XV • v 


XV • v = 


2jJ 

h4 


2y /3v v 3u u\ 

" 3h 2 \3£ + £ 3 n + nj 

^h 2 - un + vcj + xv - v 


■33 


_ 2y/v u' 


) 


^ , _ + — 1 + XV • V 

h 2 V5 n, 


T 23 ~ 0 


t 13 “ 0 


12 


h2 (| * Is 


- 2(£u + f|v) 








(A4e) 


Explicit differentiation of product terms involving the metric has been 

avoided through differentiation by parts, that is, 

^ h 3 3<J) . 8h 3 

3w h 3 3 oj + 9 3w 


For any function cj) where h^ = h 3 (w) , the term Sh-^/Sco is evaluated 
analytically. 

Equations (A4) are singular along the lines E, = 0 and r| = 0 and have 
corresponding velocities v = 0 and u = 0. To obtain the equations valid 
along these two lines, it is necessary to differentiate equations (A4) with 
respect to each of the singular coordinates, and then take the limits, 

E, -* 0 and r) -* 0, as appropriate. The resulting equations are listed below: 

Limiting form of equation along E, = 0 . - 
Continuity equation - 


3r| 


(pu) 


pu 

n 
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5-momentum equation - 


v = 0 


r|-momentum equation - 


9 ^ ( h2p u ) 


-2 IpCpuv) - |-(pu 2 + h 2 P) + (2h 2 P + pu 2 ) -2- - 


9S 


3ri 


- pu 


+ 


3t 


22 


3x 


3ri 


+ 2 




+ q {T 22 _ t 33 ) 


+ n (t 22 - t 1± ) H 


Energy equation - 


9t 


(h 2 E) = -2 ||r[v(E + P)] - |^[u(E + P)] - (E + P )(^j 


_ 3 / 3t\ 3 ( 3t\ k 3t _ 3 . . 

2 3 5\ K 35/ + 3n \ 3n) + n 3ri + 2 35 vTl1 12 


+ 3q (vT 12 + uT 22 ) + n (vT 12 + uT 22 ) 


Stress terms - 


AV • V = 


= Z2l± 


3h z 


( 2 + ** + u\ 

V 2 35 3n n/ 


: u = Sf( h2 If ~ ?v + nu) + xv - v 


t 22 = ^(» 2 ^ - nu + vC) + AV • V 


0 


33 


t 23 “ 0 


t 13 “ 0 


-12 


^ 2 (t + £0 - ^ ] 


"A 


- + if) - 2 ^] * » ' 5 


(A5b) 


) 

3 ^ t 12 

(A5c) 


(A5d) 


(A5e) 
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Limiting form of equations along p = 0.- 
Continuity equation - 


3 ,2 \ 3 . . 3 fv\ 

3^ (h f> ’ - 3E ( P V > - 2 57f'P"> - P(f) 

^-momentum equation - 

— (h 2 pv) = - g-p(pv 2 + h 2 P) - 2 ^(puv) + (2h 2 P + pv 2 )-^- - pv 

+ ?(t 1x - t 22 ) 


3t 


+ h" 


11 




+ 2 


3t 21 u_ 1 , 

^r + i (T 


ii 


t 33 ) 


ri-momentum equation - 
u = 0 

Energy equation -- 

^ - 2 U“«- + p >] 

■ <E + p, (f) + !e( k If) 

3 / 3t\ 3 . 

+ 2 3?r(f 3frj + 3C (VT 11 + UT 12 } 

+ 1‘( vT n + uT 12^ + 2 37j"( vT 12 + uX 22 ^ 


(A6a) 


) 

3 nt 2 i 

(A6b) 


(A6c) 


(A6d) 
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Stress terms - 


w ' v “ - Siltt + f + 2 


2 y ( 3v v 3u\ 

(ac + C + 9nJ 


'N 


r = 2y 
11 h 4 


^h 2 - Cv + nuj + XV • V 


22 "h 4 


2 ^ ^h 2 "577 - Du + ?v j + XV • V 


3n 


') 


2y 


• 33 " h 4 L 


- f ) - 2 - 


+ AV • V 


t 23 - 0 


T 13 “ 0 


■12 - 


h 4 t- 


,2 /^u 3v ' 

( a ? + 3n- 


2riv 


The previous dimensional equations are nondimensionalized as follows: 


P = 


_ P 


U L 


1/2 


v 


; L 


1/2 


P = 


P u 

OO CO 


e = 


A A 

C ,,T 

V oo 


€ = 


' 1/2 


n = 


>1/2 


y = 


K = 


K 

A 

Koo 


t = 


The components of the vector expression 


3w 3F(W) 3G(W) . 

3t = " 3l 3rj + H(W) + S(W) 


(A6e) 
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can now be written as follows : 


"A 



1 

CL 

(N 

1 


pv 

-> 

h 2 pv 


pv 2 + h 2 P 

w = 

h 2 pu 

F = 

puv 


h 2 E 


v (E + C 2 P) 


G = 


pu 

puv 

pu 2 + h 2 P 
u (E + C 2 P) 


H = 


-4 ♦ I) 

Q>h 2 P + p(u 2 + v 2 )] ^ - pv^- + 


[2h 2 P + p(u 2 + v 2 )]^- " P u (“ + | ) 

- (E + c 2 P)(a + |) 

where E = p [js + C]_ (u 2 + v 2 )J and 


S = 


N Re„ 


h 

5h 

h 

5n 


' h (ir + In T21 ) + 


9t> 


( E, 2 + h 2 ^ , (h 2 + 2b 2 ' 


I t 21 - f t 22 - t 33 


12 UL 22) , ( 2E, 2 + h 2 ^ , (r\ 2 + h 2 V h n 

TT ~drT/ \ fh / T 12 + V nh / T 22 n T 33 h T 11 


C 2 Y 


N. 


"a , 

L 9t\ 

9 

/ 9T \ 

^ V (1 9t 1 9t\ 

_95 

l 9?J 

+ 9p 

V an/ 

+ K U 95 n 9nJ 


Pr L 


+ ^(vt i;l + ut 12 ) + n (vT 12 + ut 22 ) 


8 9 

+ 3|‘ (vT 11 + uT 12 ) + a^ (vT 12 + UT 


22 J 


(A7) 
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n 


The W, F, and G 
0. However, the H 


terms remain unchanged along the lines 5 
and S terms are redefined as follows: 


3 (pv) _ pu 
3? n 


0 

(2h 2 P + pu 2 )-^— - pu(i) - 
- |^[v(E + C 2 P)] - (E + C 2 P)(^) 


3 (pu) _ pv 

3n 5 


H 


(2h^P + Pv 2 )-^ - pv r - 


(?) 


3n 


(puv) 


0 


for n = 0 


’(?) 


3_ 

3n 


u (E + C 2 P) 


0 

0 


3x 


22 


3x 


3n 


+ 2 


s = 


N. 


Re„ 


35 + n (T 22 T33) . 


c 2 Y 


N- 


Pr„ L 


+ 13 < T 22 " T ll ) + 3 ? T 12 

%) * U* %) + S £ 


+ 2 g^CvT-L! + ut 12 ) 


+ |^ (VT 12 + uT 22 > + ^ (VT 12 + uT 22 > 


0 and 


for 5=0 


29 


APPENDIX A 



30 


I! 



APPENDIX B 


TRANSFORMATION OF INVISCID EQUATIONS 

The conservative, vector form of the Euler equations, written in orthogonal 
parabolic coordinates (£,r),t) , was given by equation (12) as 

W t = -F jr (W) - G n (W) + H(W) 

— >• -)• ->■ _ _ t , 

where W, F, G, and H are vector quantities defined in appendix A. By using 

the transformation equations (3) 


5 = 5 


n s (£,t) - n c 

t = t 


equation (12) can be transformed from parabolic coordinates (5,Tl,t) to a set of 
coordinates in computational space (^,rj,t). To carry out this transformation, 
the following transformation operators are defined: 


9 

- 15 

3 

+ lH 


3t 

3 

35 ' 

" 35 

95 

35 

3n 

95 

3t 

9__ _ 

. ll 

9 

+ ^5- 

3 



9t 

9 

9r| 

9n 

35 

3n 

9n 

9rj 

9t 

9 

. 95 

9 

+ li 

3 

9t 

9 

9t = 

" 3t 

35 

9t 

3 rj 

9t 

dt 


(Bl) 


(B2) 


(B3) 


By using the transformation equations (eqs. (3)), 
written as follows: 


equations 


(Bl) to (B3) can be 


9_ = 9_ + in a_ 
3? 3? + 95 9n 


(B4) 


9_ = 9 jn 9_ 
9tl 9r) 


(B5) 
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5_ = 3n 3_ 3_ 

3t 3t 9r| 3t 

where 

In = I n \^s 

3C \n s - n c /3£ 


(B6) 


(B7) 


3r] _ 1 

3n " n s - n c 


3n | 

' n ^ 

,3n s 

3t ' 


I3t 


(B8 ) 


(B9) 


With the transformation operators defined by equations (B4) through (B6) , the 
vector form of Euler's equation (eq. (12)), can be transformed to computational 
space and written as follows: 


-* 


W t 


3n ■+ 

W- 

at n 





In t. 

3£ n 


in 

3n n 


-> 
+ H 


(BIO) 


Since only steady-state solutions are of interest, the term 3p/3t can be 
neglected since it approaches zero as the steady state is approached. Thus, 
the previous equation can be written as 


-> 



-Fp - F 7T- 1 G- + H 

K n an n 


where 3n/3£ and 3n/3n are defined by equations (B7) and (B8) , respectively. 
This is the same as equation (14) presented in the main text. 
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CALCULATION OF PROPERTIES AT SHOCK WAVE 

A method similar to that presented in reference 10 was the known pressure 
to calculate the shock velocity and other thermodynamic properties at the shock 
wave. Consider the shock wave illustrated in sketch (a) which can be expressed 
as n s = n s t) : 



constant 


Sketch (a) 

The shock-wave angle $ s can be expressed as 


3 S = 9 + ? s 


(Cl) 


where the slope 0 of the coordinate line r| = constant and the deviation £ s 
of the shock wave from the coordinate line ri = constant are given by the 
following equations : 


0 = tan 1 



£s = 



(C2) 


(C3) 


The derivative 3r| s /3£[ is computed numerically using central differences at 
each time step except at the downstream boundary where three-point backward 
differences are used. 
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Now, consider the velocities on the upstream side of the shock wave 
illustrated in sketch (b) : 



Sketch (b) 


The tangential component of velocity 


(V t ) ra is given by the equation 


(VtL = u oo cos e s 


(C4) 


and the normal component of velocity 


(V n )^ is given by the equation 


(V*, = U oo sin 6 s 


(C5 ) 


Changes across the shock wave can be related to the normal component, and 
the shock velocity illustrated by sketch (c) : 


(V n ) 


CO 



p 

o 

Po 

I 


Sketch (c) 

Conditions across the shock wave can be related by the following normal shock- 
wave equations : 


P 


CO 






(C6) 
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p „ + “„[< V n>„ - V s ] 2 = P s + P s [(V n ) s - V S ] 


(C7) 


(C8) 


or, since p = 1, 


YI„ 


oo V s] Ps [ 

' V „> 3 - V s] 


(C9) 

r -i 2 


r -.2 


[ ( Vco - V s] 

=■ P s + 

[<v= - 

(CIO) 

l r 

i 2 yi s 

1 r -i 2 


■ 2 [ (Vn) co - V s 

) = C 2 + 

03 1 
> 

1 

03 

> 

! | <M 

(Cll) 


Now, rearranging equation (CIO) the following expression can be obtained: 

[(Vn) s " Vs ] 2 


P ~ = [ (V n ) M - V s] i 1 - Ps 


[ (v n>„ - v s ] : 


(C12 ) 


Combining this equation with equation (C9) and solving for [( v n)„ - V s] ' the 
following is obtained: 



Ps 

Similarly, equation (Cll) becomes 



(C13) 


(C14) 


Equations (C13) and (C14) can be combined with the ideal gas equation of 


state P 



to yield the following result: 



(C15) 
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If the pressure on the downstream side of the shock wave is known, the density 
can be calculated from equation (C15) and the shock velocity can be calculated 
from the following equation: 


V 


s 



(C16) 


The velocity components in the parabolic coordinate system can now be related 
to the velocity components normal and tangent to the shock wave by the following 
equations. (See sketch (d) . ) 


U = - (V ) cos + (V t ) sin £ (C17 ) 

11 s ^ oo => 

1/ = (V n ) s sin ? s + ( V t ) m cos £ s (CIS) 



constant 
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Figure 1.- Physical and computational planes for viscous solution. 
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Figure 6.- Temperature distributions through shock layer for N Re 

with M = 5.8725 and T../T = 5.16. ' 
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Figure 7 . - Temperature distributions through shock layer for 

with = 9.92 and T /T = 10.36. 
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Figure 8.- Comparison of experimental and computed 
forebody drag coefficients C D . 


O Unpublished experimental data (CFHT) 
B Present method 
O Graves (ref. 12) 
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Figure 11.- Comparison of wall heat transfer determined 
by several different computational techniques. 
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Figure 12.- Comparison of computed wall heat transfer 
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